Annihilation dynamics during spiral defect chaos revealed by particle models

Pair-annihilation events are ubiquitous in a variety of spatially extended systems and are often studied using computationally expensive simulations. Here we develop an approach in which we simulate the pair-annihilation of spiral wave tips in cardiac models using a computationally efficient particle model. Spiral wave tips are represented as particles with dynamics governed by diffusive behavior and short-ranged attraction. The parameters for diffusion and attraction are obtained by comparing particle motion to the trajectories of spiral wave tips in cardiac models during spiral defect chaos. The particle model reproduces the annihilation rates of the cardiac models and can determine the statistics of spiral wave dynamics, including its mean termination time. We show that increasing the attraction coefficient sharply decreases the mean termination time, making it a possible target for pharmaceutical intervention. Many physical systems exhibit annihilation events during which pairs of objects collide and are removed from the system. These events occur in a number of soft-matter and active-matter systems that exhibit spatiotemporal patterning. For example, topological defects in nematic liquid crystals can develop motile topological defects that annihilate when they meet 1,2. Pair-wise annihilation of defects or singularities also plays a role in a number of biological systems. In bacterial biofilms, for instance, imperfect cell alignment results in point-like defects that carry half-integer topological charge and can annihilate in pairs. These topological defects explain the formation of layers and have been proposed as a model for the buckling of biofilms in colonies of nematically ordered cells3,4.

Pair-annihilation events are ubiquitous in a variety of spatially extended systems and are often studied using computationally expensive simulations.Here we develop an approach in which we simulate the pair-annihilation of spiral wave tips in cardiac models using a computationally efficient particle model.Spiral wave tips are represented as particles with dynamics governed by diffusive behavior and short-ranged attraction.The parameters for diffusion and attraction are obtained by comparing particle motion to the trajectories of spiral wave tips in cardiac models during spiral defect chaos.The particle model reproduces the annihilation rates of the cardiac models and can determine the statistics of spiral wave dynamics, including its mean termination time.We show that increasing the attraction coefficient sharply decreases the mean termination time, making it a possible target for pharmaceutical intervention.
Many physical systems exhibit annihilation events during which pairs of objects collide and are removed from the system.These events occur in a number of soft-matter and active-matter systems that exhibit spatiotemporal patterning.For example, topological defects in nematic liquid crystals can develop motile topological defects that annihilate when they meet 1,2 .Pair-wise annihilation of defects or singularities also plays a role in a number of biological systems.In bacterial biofilms, for instance, imperfect cell alignment results in point-like defects that carry half-integer topological charge and can annihilate in pairs.These topological defects explain the formation of layers and have been proposed as a model for the buckling of biofilms in colonies of nematically ordered cells 3,4 .

I. INTRODUCTION
In this study, we focus on the pair-wise annihilation that occurs during spiral defect chaos in excitable systems.During spiral defect chaos, spiral waves continuously break down to form new ones and are removed through collisions with other spiral waves.While spiral defect chaos is present in a variety of chemical and biological pattern-forming systems [5][6][7][8][9][10][11][12][13][14][15][16] , perhaps its most studied example can be found in models of cardiac tissue [17][18][19][20] .These models naturally exhibit spiral waves and the tips of these spiral waves undergo stochastic annihilation events [21][22][23] .Importantly, these annihilations underlie cardiac fibrillation, characterized by unorganized electrical wave propagation in the heart 20 .Fibrillation in the ventricles is lethal 24 while atrial fibrillation, the most common cardiac arrhythmia in the world with approximately 30 million patients in 2010, is associated with increased morbidity and mortality [25][26][27][28] .
Previous studies have described the statistics of spiral wave tips in spatially extended systems in the form of a stochastic birth-death process 11,13,29 .More recently, this approach was applied to cardiac models, where the creation and annihilation rates of spiral tips were determined for various domain sizes 22,30 .Since N = 0 is an absorbing state, the dynamical state will eventually terminate.Using these rates, the mean termination time, τ, was computed and was shown to be exponentially distributed, consistent with experiments and clinical data 22,30,31 .This termination time is a quantity of interest in the context of cardiac dynamics as termination indicates the heart has transitioned into normal sinus rhythm.Thus, minimizing τ is of critical importance for managing cardiac fibrillation.
Previous work has shown τ depends on the tissue size, A, and reducing τ can be achieved by decreasing A 14,22 .Unfortunately, decreasing the effective size of cardiac tissue is not practical and determining the dependence of τ on other physiological properties is therefore desirable, especially if these properties can be drug-targeted.In this study, we propose targeting the attraction coefficient, a, which controls the strength of attraction between spiral wave tips.
We first simulate two spatially extended cardiac models and show that the spiral tips that annihilate in these models display an apparent attractive interaction at short distances and diffusive Brownian behavior at large distances.We then formulate a stochastic particle model in which tips are represented as particles and show that it can capture the attractive and diffusive properties of the tips in the cardiac models.Furthermore, we show that this particle model generates tip dynamics that reproduce both the annihilation rates as a function of the density of tips and the distribution of termination times for the two cardiac models.Finally, we show that increasing the attraction coefficient significantly decreases the mean termination time of spiral defect chaos.

II. CARDIAC MODELS
To determine the dynamics of spiral wave tips in the cardiac models, we integrated the mono-domain equations, which describe the time evolution of the membrane potential, u, by the excitable reaction-diffusion equation where I ion represents the transmembrane currents, C m = 1µF/cm 2 is the transmembrane capacitance per unit area, and D u is the diffusion coefficient 20 .To stress the generality of our approach, we used two commonly-employed models for cardiac tissue: the Fenton-Karma (FK) model 32 and the Luo-Rudy (LR) model 33 .The precise formulation of I ion for these models is provided in Appendix A, along with model parameters.
We integrated Eq. 1 on a square computational domain of size A and enforced periodic boundary conditions.These periodic boundary conditions resulted only in pair-creation and pair-annihilation of spiral tips due to conservation of vorticity in u and, thus, in even numbers of tips.We used a spatial discretization of ∆x = 0.025 cm and a temporal discretization of ∆t = 0.025ms and used a body-centered forward-time explicit Euler method with the Laplacian operator discretized using a five-point stencil on the square lattice.
We chose parameter values for which both models exhibit spiral defect chaos.In other words, a single spiral wave is unstable and will break up into multiple spiral waves.Observations of spiral tip motion began 100 ms after the start of the simulation at time t = 0 so as to allow periodic boundaries enough time to become smooth.The locations of spiral tips were determined from the intersection points of two level sets of constant voltage (u = 0.4 for the FK model and u = −30 mV for the LR model).One level set was the equipotential line at the constant threshold voltage.The other level set was the equipotential line at the constant threshold voltage 4 ms later in time.Linear segments of level sets were determined by linear interpolation along the edges that connected nodes using the marching squares method described by Lewiner et al. 34 modified to support the periodic boundary conditions.Intersection points were determined by solving the linear system of equations that describes two intersecting linear segments and were recorded as spiral tip locations.For the chosen parameter values, the width of a planar wave was 0.41 cm for the FK model and 0.15 cm for the LR model, computed at the same thresholds used for the detection of spiral tip locations.

III. CARDIAC MODEL RESULTS
Snapshots of typical simulations displaying a pair-wise annihilation event are shown in Fig. 1A for the LR model (upper row) and the FK model (lower row).In these panels, u is visualized using a grayscale and clockwise and counterclockwise tips are indicated by black and yellow symbols, respectively.Examples of single spiral tip trajectories are shown in Fig. 1B, with the FK trajectories plotted in blue and the LR trajectories in orange.Although spiral tips can be both created and annihilated, we focus here on the behavior of annihilating spiral wave tips.In our simulations, only pairs of counter-rotating spiral waves that are connected by an activation front can annihilate.Annihilation occurs when a depolarized region acts as a wave block, causing the activation front to shrink in length before spontaneously disappearing.
For both models, we examined whether tips that came within a certain distance were able to escape and move apart or always annihilated.Specifically, we determined the number of pairs that became separated by less than 0.1 cm and were then able to move further apart than 0.15 cm.These pairs always consisted of spiral waves of opposite chirality: one was rotating clockwise while the other was rotating counterclockwise 35,36 .Examining 355 pairs in the FK model and 99 pairs in the LR revealed that these pairs always annihilated.In the remainder of this study, we will strictly focus on annihilation events and all subsequent results are therefore conditioned on these dynamics.I).
To quantify the dynamics of annihilating spiral tips, we first computed the lifetime of pairs that annihilate, Γ.For this, pair-annihilation events were determined by ordinary nearestneighbor particle tracking subject to periodic boundary conditions.We computed Γ as the time of annihilation minus the time of creation of the younger of the pair.The distribution of these lifetimes was approximately exponential (Fig. 2A) with the FK model producing significantly longer-lived pairs on average (Table I).
The distribution of termination times of spiral defect chaos, i.e. the time until all spiral tips have been removed, was also found to be fitted well with an exponential (Fig. 2B, see also 22 ).This can be understood by considering not only pair annihilation but also pair creation, which may occur when a wave back meets a wave front.As a result of this contin- uous creation and annihilation process, the number of tips, N, reaches a long-lived metastable state 37,38 .The distribution associated with this metastable state is called the quasistationary distribution and can be calculated using the master equation for the probability P(N,t) to have N = 0, 2, 4, 8, ... tips at time t 22 : where W ±2 are the annihilation and creation rates.If the number of spiral tips reaches zero, the spiral defect chaos has terminated.In other words, N = 0 in the above equation is an absorbing boundary, and W +2 (0) = 0. Through repeated simulations of termination events, we determined the distribution of termination times, which is shown for both models in Fig. 2B.By averaging 6,043 tip trajectories at a domain size of A = 25 cm 2 , we next computed the mean squared displacement (MSD) as a function of time lag 39 (Fig. 3A).For both cardiac models, the MSD was not significantly different from linear for long timescales with exponent values of 1.002±0.012(FK model for time lags from 100ms to 4000ms) and 0.996±0.019(LR model for time lags from 20 ms to 200 ms).Multiple other fits using different time lag windows supported this finding (Fig. 4).
Thus, in both cardiac models, the spiral wave tips can be described as effectively undergoing Brownian motion for long timescales.At short timescales, however, tips did not exhibit diffusive behavior.To determine the behavior of spiral tips at these short timescales, we analyzed the movement of annihilating spiral waves in the simulations.This analysis revealed that annihilation occurs when the activation front connecting the tip pair is blocked by a nearby polarized region.The block results in a rapid shrinking of the activation front and the removal of the pair (see Movie S1).Thus, at short timescales, the annihilating pair of tips appear to attract, which becomes apparent from Fig. 3B where we plot the ensemble-averaged radial velocity dR/dt as a function of the distance between the tips, R.This velocity is not constant, but instead becomes more negative when 1/R increases (and R decreases).In other words, an apparent attractive force induced the annihilation of counter-rotating pairs of spiral waves.We should stress that this force is only apparent and does not arise from a physical attraction between the different spiral waves.
We next sought to quantify this attractive force.Since we were interested in the attractive force during annihilation, we only focused on pair-wise annihilation events in our simulations, ignoring creation events.We identified 51,452 annihilating pairs and, for a given pair annihilating at time t f , we determined R versus time until annihilation, t ′ ≡ t f − t ≥ 0. From this, we computed the mean squared range (MSR) by ensemble-averaging R 2 conditioned on a given t ′ .The results of this ensemble-averaging are shown in Fig. 5 as solid lines, together with the 95% confidence intervals as shaded areas.On short timescales, the MSR demonstrates oscillations, illustrating that the apparent attractive force has an oscillatory component.As will be discussed, this oscillatory component is simply the manifestation of the rotation of the spirals.FIG. 5. MSR between annihilating tips versus time until annihilation from simulations of the FK and LR models with shaded regions corresponding to 95% confidence intervals.Also shown are the fits of MSR from the OPM (solid lines) and the LPM (dashed lines).

IV. STOCHASTIC PARTICLE MODELS A. Oscillatory Particle Model
To gain further insight into the annihilation dynamics of the tips using the spatially extended cardiac models is computationally expensive, especially for large domain sizes.Therefore, we developed a model in which the spiral wave tips are represented by moving particles subject to an oscillatory short-range attractive force and Brownian motion with diffusion coefficient, D. Such a significant simplification of a cardiac model was earlier used to describe the chaotic tip trajectories of a single spiral wave in the presence of heterogeneities 40 .
In this oscillatory particle model (OPM), we modeled the attractive force between two annihilating tips as inversely proportional to the distance between them.This attractive force was assumed to consist of a constant term, parameterized by the coefficient a 0 , together with an oscillatory term, parameterized by a 1 .This oscillatory term arises from the fact that each spiral wave rotates and that the attractive force depends on the relative location of the two activation fronts 35,36 .This relative location is determined by the angle of the activation front (AAF), θ , with an arbitrary, fixed axis of each spiral wave.For example, one could choose to measure θ relative to the positive x-axis.Measuring this AAF, however, is challenging since it requires determining the exact location of both the tip and the activation front and is thus prone to large discretization errors.For this reason, we will not explicitly quantify the AAF in this study.
Based on observations of the cardiac model simulations, we will take the frequency ω and period T OPM of the annihilating pair to be identical.Furthermore, as we pointed out before, the spiral waves of an annihilating pair have opposite chirality so that AAF of the first particle progresses as θ 1 (t) = ωt + θ 1 and that of the second one as θ 2 (t) = −ωt + θ 2 .Then, the difference between the two AAFs can be written as ∆θ /2 = ωt + θ f , where θ f = (θ 1 − θ 2 )/2.This difference between the two AAFs is maximal and minimal once during a complete rotation of the two tips.Based on the oscillatory nature of the MSR curves in Fig. 5, we will assume that the attractive force depends on this difference.Taken together, the deterministic part of the OPM for one of the particles, located at (X 1 ,Y 1 ), is written as 2 is the distance between this particle and the other particle, located at (X 2 ,Y 2 ).The deterministic OPM equations for this other particle are then given by dX 2 dt = − dX 1 dt , dY 2 dt = − dY 1 dt .To capture the diffusive motion of the spiral tips, we included an independent Brownian motion term to the deterministic equations for each particle.Then, we find that the distance between annihilating tips can be modeled by the following Langevin equation: where W t is a Wiener process at time t and D is estimated from the least-square slope of MSR at large ranges (MSR > 3 cm 2 ).Note that the factor of two arises from the pair-wise interaction.Furthermore, since the variance of the two twodimensional tips add, the Wiener process is multiplied by √ 8D instead of √ 2D.Using the Langevin equation, it is straightforward to derive an expression for the MSR, < R 2 >.This MSR has both a deterministic and a stochastic term: This result was verified by comparing this expression to the average MSR of 10,000 statistically independent explicit simulations of Eq. 5.These simulations can be carried out by either starting with a large value of R and progressing until R falls below a small threshold value or by starting with a small value for R and integrating backwards in time.The next step was to fit Eq. 6 to the MSR curves obtained from the cardiac models using simulated annealing on the last 300 ms (FK) and 100 ms (LR) before annihilation.The fits to both cardiac models for these time intervals are excellent (mean percent error (MPE) <4%), as can be seen Fig. 5.The fitted parameters, a 0 , a 1 , θ f , and T OPM are reported in Table I, together with the aforementioned estimates for D. In Table I we also report the period of the spiral waves in the cardiac models, T , determined by computing the mean number of rotations per lifetime.A comparison between T OPM and T reveals that these periods match perfectly, indicating that the oscillatory component in the MSR is due to the rotational nature of the spiral wave.

B. Linear Particle Model
The OPM can, in principle, be used to compute annihilation rates, which can then be compared to results from the cardiac models.However, computing these rates with the OPM is challenging because annihilation only occurs when the final difference in AAF between the two particles takes a specific value.Thus, it requires tediously tracking the AAF for each particle and enforcing the final condition θ (t f ) = θ f at the time of annihilation.To circumvent this problem, we simplified the OPM to a linear particle model (LPM) by taking the attractive force between particles as linear in 1/R and dropping the oscillatory term.In this case, the attractive force is parameterized by a single coefficient a and, modeling the diffusive behavior of spiral tips as before, the LPM reads which results in a MSR given by To relate a to the parameters of the OPM, we demanded that the MSR averaged over the exponentially distributed lifetimes, ⟨MSR⟩, be equal for both particle models.For the LPM, we can derive while a similar calculation for the OPM gives (10) Setting ⟨MSR LPM ⟩ = ⟨MSR OPM ⟩ results in an analytical expression for a in terms of the OPM parameters that is independent of D: The estimates of a evaluated from Eq. 11 are listed in Table I, and corresponding MSR LPM plots are shown as dashed lines in Fig. 5. Repeating the analysis for different domain sizes, A, revealed that the estimate of a was largely independent of A for both of the cardiac models.An example of the MSR obtained using A = 39.0625cm 2 is shown in Fig. 6A.For this, and all other domain sizes, we were able to fit the MSR using the oscillatory particle model (OPM), as demonstrated by the solid lines in Fig. 6A, which allowed us to compute the value of a.The resulting value of a decreased slightly as the domain size increased (Fig. 6B), suggesting that the LPM can be used to simulate spiral tip annihilation at different domain sizes using a fixed set of model parameters determined for a single value of A. Furthermore, the sum a + 2D also changed little as the domain size increased (Fig. 6C).The product of this sum and the mean lifetime is proportional to the mean squared distance (see Eq. 9).These results indicate that the average distance between annihilating particles is determined by local properties, which are largely insensitive to the domain size.
To compute annihilation rates in using particle dynamics instead of cardiac model simulations, we implemented the LPM, using the obtained values of a and D, by simulating N particles with locations (X 1 (t),Y 1 (t), . . ., X N (t),Y N (t)) on a square domain of size L 2 with periodic boundary conditions.Each particle i is time stepped with time step ∆t = 10 −2 ms and the X coordinate of particle i obeys where D > 0 is the diffusion coefficient, Z is a random value with zero mean and unit variance, and R i j is the distance between particle i and j.A similar equation can be written for Y i .
For initial conditions, we considered N uniformly distributed particles at two domain sizes (A = 25 cm 2 and A = 100 cm 2 ) and updated particle positions every ∆t = 0.01 ms.Pairs of particles were removed from the simulations, and thus annihilated, at rate κ whenever they were closer than a reaction range, r.As an estimate for this reaction range, we chose the 25 th percentile of the distribution of closest distances between non-annihilating tips.We have verified, however, that other choices of r give similar qualitative results.

A. Annihilation Rates
We used the LPM to compute annihilation rates and adjusted the only free parameter (κ) to match the annihilation rates found in the cardiac models.This is facilitated by the fact that in our previous work we showed that the latter, computed for different values of A, collapsed onto a single curve when rescaled by A 22 .Our results show that the annihilation rate can be approximated by a power law, w − = M − n ν − , where w − ≡ W − (N)/A is the rescaled annihilation rate and n = N/A is the tip density (symbols Fig. 7).The fitted LPM annihilation rates are shown as dashed lines in Fig. 7 and resulting parameters are listed in Table I.These fits demonstrate that the LPM can accurately replicate the annihilation rates of the cardiac models (MPE <4%).Importantly, as was the case for a, these fits use the same parameter values for both system sizes.Also note that simulations of the particle model are much more efficient than the cardiac models, especially for large domain sizes.Specifically, the particle model simulations use O((L/∆x) 2 − N 2 avg ) fewer operations per time step, where N avg is the average number of particles.For example, for A = 100 cm 2 , the speed-up exceeded 10 4 -fold.FIG. 7. Mean annihilation rate versus number density for spiral tips from the cardiac models (symbols) and their linear particle model fits (dashed lines).

B. Varying Attraction Strength
Although other parameters of the LPM can be varied as well, we focused on the sensitivity of the annihilation rate to the model parameter a, which controls the strength of attraction between particles.By carrying out additional simulations of the LPM, we determined how the annihilation rate depends on a while holding the remaining parameters D, r, and κ fixed.This rate was found to increase with increasing values of a, which can be understood by realizing that larger attractive forces result in distant particles coming closer together faster.Importantly, however, we found that the rate was always fitted well by the power law w − = M − n ν − (Fig. 8A).For both models, we found that for increasing values of a the exponent ν − became smaller (Fig. 8B) while M − increased (Fig. 8C).

C. Statistical Properties
Using the fact that the annihilation rate in the LPM can be fitted by a power law, we were able to compute statistical properties of spiral defect chaos in the spatially extended car-diac models.For this, we used the creation rate in the cardiac models, which were computed by counting the rate of creation events conditioned on number density.As with the annihilation rates, we found the creation rates to be captured by a power law fit, as well: w + = M + n ν + (Fig. 9A).In the remainder, we will keep this creation rate fixed while varying a, and thus the annihilation rate, according to Figs. 8B&C.Using the power laws for creation and annihilation rates, we constructed a transition matrix from which we computed the distributions of termination times 22,41 .These distributions were found to be exponentially distributed for all values of a (Fig. 9B), consistent with experimental data 31 .Furthermore, the distribution tends to smaller termination time values for increasing a.It is now also possible to obtain approximate expressions for quantities relevant to cardiac models and fibrillation in clinical settings.For example, the average number of particles, N avg can be estimated by solving w + (n) = w − (n), resulting in Thus, the mean particle number is predicted to increase linearly with the domain size.Fig. 10A shows the results from Eq. 12, using parameters corresponding to A = 25 cm 2 for all domain sizes, which is possible since a is insensitive to A (Fig. 6B).The estimated value using Eq. 12 agrees well with the average number of tips from the cardiac models (symbols), especially for larger values of A. Furthermore, the average number of tips computed using Eq. 12 decreases as a function of a (Fig. 10B), which can be explained by the fact that ν − decreases for increasing values of a (Fig. 8B).
In addition, the mean termination time can be computed using an analytic solution 42 .Specifically, for paired birth-death processes, the mean termination time τ is a function of the initial (even) number of spiral tips, N 0 : where ϕ(2k) = ∏ k i=1 w − (2i/A)/w + (2i/A) and ϕ(0) ≡ 1.Since the termination time is dominated by τ(2), we can use the power laws for the annihilation and creation rates to derive an approximate closed expression.For large values of A, the termination time can be shown to be proportional to 22 which becomes, after substituting the expressions for w ± , Thus, consistent with earlier work using the cardiac models 14,22 , the termination time increases exponentially with A. Importantly, aside from a prefactor that is independent of A, this expression provides an explicit estimate for τ for any domain size, including ones that would be prohibitively expensive to simulate directly.In Fig. 10C, we plot the values for τ obtained from the cardiac models as symbols together with a fit using Eq.14.The former ones were, of course, only obtainable for small domain sizes while the fit can be extended to arbitrarily large domain sizes.
Finally, using the fitted power laws (Fig. 9B), we computed τ as a function of a with the other parameters fixed to those corresponding to the cardiac models.We found that τ decreased by a factor of ∼ 10 2 in response to a increasing by a mere factor of 10, as shown in Fig. 10D.This great sensitivity of τ to a, considered together with the relative insensitivity of a to A, suggests that changing a is a mechanism for controlling the mean termination time for minimally varying domain sizes.Specifically, if modifying the electrophysiological parameters of the cardiac models increases attraction only, then τ decreases, resulting in shorter termination times.
To verify that it is possible to vary a, and thus τ, we carried out additional simulations of the cardiac models.For this, we altered the excitability parameter τ d in the FK model and the extracellular potassium concentration [K + ] o in the LR model.We followed the same procedure as before and fitted the MSR expression for the OPM to the MSR curves obtained from the cardiac models using simulated annealing on the last 300 ms (FK) and 100 ms (LR) before annihilation.The fit results are shown in Figure 11 for both models.These fits, together with measured values of D and Γ, were then used to compute a and the LPM model was simulated to determine τ.
For the FK model, we increased τ d by 20% and obtained as fitted parameter values a 0 = 1.6077 cm 2 /s, a 1 = 1.3311 cm 2 /s, T OPM = 110.3665ms, and θ f = −1.0730radians (RMSE=0.0459cm 2 ).The measured parameters were D = 0.048 ± 0.007 cm 2 /s and Γ = 101.3±1.3 ms, resulting in a = 1.82 ± 0.01 cm 2 /s, which was 17.5% larger than the corresponding value for a in Table I.This led to a decrease of the mean termination time by 35%.In the LR model, we reduced [K + ] o by 29.6% and found the following fitted parameter values: a 0 = 4.1932 cm 2 /s, a 1 = 5.9593 cm 2 /s, T OPM = 82.4087ms, and θ f = −0.9490radians (RMSE=0.0352cm 2 ).Measured parameters were D = 0.42 ± 0.03 cm 2 /s and Γ = 45.9 ± 1.1 ms, resulting in a = 5.73 ± 0.07 cm 2 /s, which was 38.6% smaller than the corresponding value for a in Table I.The termination time in this case increased by 173.6%.Thus, increasing a reduced τ for the FK model, and reducing a increased τ for the LR model.

VI. CONCLUDING REMARKS
In summary, this study reveals that the annihilation dynamics of spiral waves in spatially extended cardiac models can be captured by a computationally efficient particle model.In this model, the spiral wave tips are represented by diffusing particles subject to a locally attractive force.We showed that the particle model accurately reproduced the annihilation rates of spiral tips in the cardiac models and that it can be used to efficiently model their tip dynamics.We also used the particle model to show that increasing the strength of the apparent attraction force accelerates annihilation, thus decreasing the mean termination time.
What have we ignored and how can this study be extended?First of all, we exclusively focused on annihilation dynamics.Clearly, to construct more realistic particle models, we will have to include creation events.These events may require the inclusion of repulsion of spiral wave tips, as has been reported in some cardiac models 35 .Once these events are incorporated into the model, we should be able to determine how parameters that describe these creation events can change termination times.Second, we have simulated the dy-namics on a domain with periodic boundary conditions, guaranteeing that annihilation is a pair-wise event.Introducing non-conducting boundaries will necessitate a description of the annihilation of a single spiral wave that collides with one of these boundaries.Third, we have simplified the oscillatory model to a linear model such that we did not have to compute the AAFs for the spiral waves.It would be interesting to extend our study to explicitly enforce conditions on the relative AAFs at the time of annihilation.Finally, future work involving cardiac models can include a more thorough investigation of the dependence of the attraction coefficient between spiral tips on trans-membrane currents.Reproducing desirable effects on these currents can then be the target of drug discovery, potentially opening a new door to noninvasive therapies for clinically significant symptoms of cardiac fibrillation.

VII. ACKNOWLEDGEMENTS
This work was supported by the National Institute of Health through awards HL083359, HL122384, HL149134, and T32GM127235.T.T. acknowledges the Open Science Grid 43,44 for allocation of computational resources.P.M acknowledges the support of the Research Training in Mathematical and Computational Biology Grant (1148230) and the NSF Division of Mathematical Sciences award 1903275.The ionic currents in the Fenton-Karma (FK) cardiac model 45 are written as: where the fast inward current is the slow outward current is and the slow inward current is In these equations, v and w represent the fast and slow variable, respectively: Parameter values for the model were chosen to be τ d = 0.45ms, τ 0 = 12.5ms, τ r =33.25ms, τ si = 29ms, k = 10, u si c = 0.85, u c = 0.13, u v = 0.04, τ − v1 = 1250ms, τ − v2 = 19.6ms,τ + v = 13.03ms,τ − w = 40ms, and τ + w = 800ms.Using a diffusion coefficient of D u = 0.0005cm 2 /ms resulted in a conduction velocity of a planar wave of c v = 51cm/s, which is within the electrophysiological range.

Luo-Rudy Model
The total current for the Luo-Rudy (LR) cardiac model 46 used in this study is given by where the first term is the sodium current with parameter values G Na = 16mS/cm 2 and E Na = 54.4mV.The second term is the slow inward current where G si = 0.052mS/cm 2 and E si = −82.3mV− 13.0287 log [Ca +2 ] i ) .The last two terms represent the potassium current, which is computed using the equations and parameters of Qu et al. in Ref. 47,48 .In these expression, the intracellular calcium is given by In these equations, a y , b y are dimensionless monotonic functions, which were evaluated in constant time using a lookup table based on the equations of Qu et al. 47,48 .Additional parameters can be found in these references.For complete reproducibility, the equations and parameters that we used are included in , and where we have taken [K + ] o = 5.4mM.We computed I K using I K = I 1 x, where where have used the values G K = 0.423mS/cm 2 and where we took [K + ] i = 145mM, [Na + ] i = 18mM, [Na + ] o = 145mM, R = 8.3145J/(mol • K) as the universal gas constant, and F = 96485.3321233100184C/molas Faraday's constant.We supposed a homeostatic body temperature of T = 37 • C fixed.We have similarly approximated To avoid numerical overflow when time evolving these gating variables, we treated time scales as zero, τ y = 0, whenever they took a value τ y ≤ 5 • 10 −4 ms.Additionally, we evaluated I 1 using L'Hospital's Rule in the 10 −6 mV neighborhood of V m = −77mV.
Finally, for the diffusion coefficient we used D u = 0.001cm 2 /ms, resulting in a conduction velocity of c v = 33cm/s, which is within the electrophysiological range.TABLE III.Parameters of the Luo-Rudy-I model used in our study.For ease of notation, we used the dimensionless transmembrane voltage field u = V m /(1mV ).Gating variables time evolve according to Eq. SA3.

FIG. 1 .
FIG. 1.A Grayscale snapshots of u showing spiral defect chaos in (top) the LR model and (bottom) the FK model with A = 25 cm 2 .Indicated are the tips of clockwise (black stars) and counterclockwise (yellow stars) rotating spiral waves.Snapshots were taken at (left) t ′ = 8 ms, (middle) t ′ = 4 ms, and (right) t ′ = 0 ms before an annihilation event.Annihilation can be explained by a wave-block resulting from a depolarized area acting as a wall to spiral tip motion.B Single tip trajectories are shown from the FK model and the LR model (orange).The trajectories in the LR model show fewer pivots, supporting our result that the LR model has a larger diffusion coefficient than the FK model (see also TableI).

FIG. 2 .
FIG. 2. A Probability density of lifetime Γ of annihilating pairs in the two cardiac models.Inset is the probability density of Γ to visualize the short-lived tips, revealing that the abundance of shortlived spiral tips (< 15ms) was somewhat greater for the LR model relative to the FK model.This log-log histogram has 10 bins per decade.B Probability density of the termination times of simulations of the two cardiac models, along with an exponential fit (solid line).
FIG. 6.A MSR between annihilating tips versus time until annihilation from simulations of the FK and LR models, using a larger domain size of A = 39.0625cm 2 , with shaded regions corresponding to 95% confidence intervals.The solid lines correspond to fits from the OPM while the dashed lines correspond to fits from the LPM.B Computed attraction coefficient versus domain size.C Sum of attractive and diffusive forces versus domain size.Error bars indicate 95% confidence.

AFig. 3 FIG. 8 .
Fig. 3 Results of particle models when changing a \textbf{a} Mean annihilation rate density versus number density.Power laws were fit from the linear particle model for increasing values of the attraction coefficient, $a$, as indicated by (inset) the color bar.\textbf{b} Power law exponent as a function of $a$ for (blue) the FK model and (orange) the LR model.\textbf{c} Power law magnitude versus $a$.

FIG. 9 .
FIG.9.A Mean creation rate (triangles) and mean annihilation rate (dots) versus number density for spiral tips from the cardiac models using different domain sizes.Dashed lines correspond to power law fits (TableII).Black dots correspond to the mean particle density.B Probability density of termination times of the LPM for increasing values of a equally spaced from a = 1 cm 2 /s to a = 5 cm 2 /s.Parameter values correspond to the FK model and A = 25 cm 2 .

DFIG. 10 .
FIG.10.A Average tip number as a function of A computed using the cardiac models (symbols), along with the linear prediction of Eq.S12.B Average tip number as a function of a computed using the LPM with parameter values corresponding to A = 25 cm 2 .The darkened symbols correspond to the value of a representing the cardiac models.C Mean termination time versus A computed using Eq.S15 (dashed lines) and separately obtained from the cardiac models (symbols).D Mean termination time as a function of a (using parameter values for A = 25 cm 2 ).Black circles correspond to a obtained by fitting the cardiac models.

FIG. 11 .
FIG. 11.A MSR from the FK model with the excitability parameter, τ d , increased by 20%.The blue shaded region represents 95% confidence intervals estimated via bootstrap.The black line represents the fit to the OPM.B MSR from the LR model with the extracellular potassium concentration, [K + ] o , decreased by 29.6%.The orange shaded region represents 95% confidence intervals estimated via bootstrap.The red line represents the fit to the OPM.

TABLE I .
Particle properties of spiral tips from the cardiac models including parameter values corresponding to OPM and LPM.

TABLE II .
Power law fits for the annihilation rates and the creation rates.Also included is the mean termination time and average number of particles.